{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 二维虚单元法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 记号\n",
    "| 记号 | 意义 |\n",
    "| -------|---------|\n",
    "|$\\mathbb{R}^2$ | 二维欧氏空间 |\n",
    "|$\\mathcal{D}$ | 二维区域 |\n",
    "|$\\bf{x}_{\\mathcal{D}}$ | 区域 $D$ 的 重心 |\n",
    "|$\\mathcal{P}_k(\\mathcal{D})$ | 区域 $D$ 上 $k$次多项式空间| \n",
    "|$E$ | 多边形区域  |\n",
    "|$N^V$ | $E$ 的顶点个数 |\n",
    "|$\\bf{x}_i$ |  第 $i$ 个顶点坐标|\n",
    "|$|E|$ | $E$ 的面积 |\n",
    "| $h_E = \\sqrt{|E|}$ | $E$ 的尺寸 |\n",
    "|$\\mathbf W=\\begin{pmatrix}0&-1\\\\1 & 0 \\end{pmatrix}$ | 旋转矩阵, 作用在二维列向量上表示逆时针旋转该向量$90^\\circ$ |\n",
    "|$\\alpha=(\\alpha_1, \\alpha_2)$ | 二重指标, $\\mathbf x^\\alpha = x^{\\alpha_1}y^{\\alpha_2}$|\n",
    "|$\\mathcal{M}_k(\\mathcal{D})$ | 单项式空间 |\n",
    "|$\\mathcal V_k(E)$ | $k$ 次的虚单元空间 |  \n",
    "|$\\Pi^\\nabla: \\mathcal V_k(E) \\rightarrow \\mathcal{P}_k(E)$ | $k$ 次虚单元空间到 $k$ 次多项式空间的的投影算子| \n",
    "|$\\Pi^\\nabla_*: \\mathcal V_k(E) \\rightarrow \\mathcal V_k(E) $ | 投影算子  | \n",
    "|$ n_k = (k+1)(k+2)/2$ | $k$ 次多项式空间维数 |\n",
    "|$ N_k = \\text{dim} \\mathcal V_k(E) = N^Vk+(k-1)k/2$ | $k$ 次虚单元空间的维数 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  模型\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "-\\Delta u = f, & \\text{ in } \\Omega \\\\\n",
    "u = g, & \\text{ on } \\partial\\Omega \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "找到 $u \\in \\mathcal H^1_g(\\Omega)$:\n",
    "$$\n",
    "(\\nabla u, \\nabla v) = (f, v),\\quad\\forall v \\in \\mathcal H^1_0(\\Omega)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 梯度算子的矩阵表达"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla (m_1, m_2, ..., m_{n_k}) = (m_1, m_2, \\cdots, m_{n_{k-1}})(G_x, G_y)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Laplace 算子的矩阵表达"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Delta (m_1, m_2, \\cdots, m_{n_k}) = (m_1, m_2, \\cdots, m_{n_{k-2}})L\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "综上可得, 当 $a < 2$,且 $b<2$ 时, 有:\n",
    "\n",
    "$$\n",
    "\\Delta(a,b) = 0\n",
    "$$\n",
    "\n",
    "当 $a \\geq 2$,或 $b\\geq 2$ 时, 有:\n",
    "\n",
    "$$\n",
    "\\Delta(a,b) = \\frac{\\max\\{a,b\\}(\\max\\{a,b\\}-1)}{h_E^2}[(a-2,b) + (a,b-2)]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里规定 $a-2<0$ 时, $(a-2,b) = 0$; $b-2<0$ 时, $(a,b-2) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "若已知 $(a,b)$, 则对应的 $\\alpha$ 编号为:\n",
    "\n",
    "$$\n",
    "\\alpha = \\frac{(a+b)(a+b+1)}{2}+b\n",
    "$$\n",
    "\n",
    "若已知 $\\alpha$, 则对应的 $(a, b)$, 计算公式为\n",
    "$$\n",
    "a+b = \\lfloor \\frac{-1+\\sqrt{1+8\\alpha}}{2}\\rfloor\\\\\n",
    "b = \\alpha - \\frac{(a+b)(a+b+1)}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 虚单元空间 $\\mathcal V_k(E)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于任一 $v_h\\in\\mathcal V_k(E)$, $v_h$ 满足如下性质：\n",
    "\n",
    "* $v_h|_e \\in \\mathcal P_k(e)$\n",
    "* $v_h|_{\\partial E} \\in \\mathcal C^0(\\partial E)$\n",
    "* $\\Delta v_h \\in \\mathcal P_{k-2}(E)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知 $\\mathcal V_k(E)$ 的维数为 \n",
    "\n",
    "$$\n",
    "N_k = \\text{dim} \\mathcal V_k(E) = N^Vk+ \\frac{(k-1)k}{2}.\n",
    "$$\n",
    "\n",
    "设其基函数集合为 $\\{\\varphi_i\\}_{i=1}^{N_k}$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$v_h\\in \\mathcal V_k(E)$ 的自由度定义如下：\n",
    "\n",
    "* $v_h$ 在 $E$ 的所有顶点处的函数值\n",
    "* $v_h$ 在 $E$ 的每条边内部 $k-1$ 个 Gauss-Lobatto 积分点处函数值\n",
    "* $v_h$ 在单元 $E$ 内部与 $\\mathcal P_{k-2}(E)$ 中所有单项式 $m_\\alpha$ 的乘积的积分平均值， $$ \\frac{1}{|E|}\\int_E v_h m_\\alpha, \\quad \\alpha = 1, \\ldots, n_{k-2}$$\n",
    "其中 $n_{k-2} = \\dim \\mathcal P_{k-2}(E)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 投影算子\n",
    "\n",
    "定义 $\\Pi^\\nabla: \\mathcal V_k(E)\\rightarrow \\mathcal{P}_k(E)$ 如下: 给定  $ v_h \\in\n",
    "\\mathcal V_k(E)$, 找到 $\\Pi^\\nabla v_h \\in \\mathcal{P}_k(E)$ 满足:\n",
    "\n",
    "$$\n",
    "(\\nabla \\Pi^\\nabla v_h, \\nabla p_k) = (\\nabla v_h, \\nabla p_k), \\forall p_k\\in\\mathcal{P}_k(E)\n",
    "$$\n",
    "\n",
    "投影 $\\Pi^\\nabla v_h$ 关于常数唯一.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以定义投影 $ P_0: \\mathcal V_k(E) \\rightarrow\n",
    "\\mathcal{P}_0(E)$:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "P_0(\\Pi^\\nabla v_h - v_h) = 0\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "具体定义如下: \n",
    "$$\n",
    "\\begin{align}\n",
    "P_0 v_h :=& \\frac{1}{N^V}\\sum_{i=1}^{N^V} v_h(\\mathbf x_i)\\text{, for } k=1\\\\\n",
    "P_0 v_h :=& \\frac{1}{|E|} \\int_{E}v_h = \\frac{1}{|E|} (1, v_h)_E\\text{, for } k\\geq 2\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "其中 $\\frac{1}{|E|} (1, v_h)_E$ 为 $v_h$ 在 $E$ 内部自由度处的值."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Pi^\\nabla \\varphi_i = \\sum_{\\alpha = 1}s_{i,\\alpha}m_\\alpha\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "(\\nabla \\sum_{\\alpha = 1}^{n_k}s_{i,\\alpha}m_\\alpha, \\nabla m_\\beta) = (\\nabla \\varphi_i, \\nabla m_\\beta), \\quad \\beta = 1, \\ldots, n_k\\\\\n",
    "\\sum_{\\alpha = 1}^{n_k}s_{i,\\alpha}(\\nabla m_\\alpha, \\nabla m_\\beta) = (\\nabla \\varphi_i, \\nabla m_\\beta), \\quad \\beta = 1, \\ldots, n_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以得到上述投影的矩阵表示形式:\n",
    "$$\n",
    "\\mathbf G \\mathbf S = \\mathbf B\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf S = \n",
    "\\begin{pmatrix}\n",
    "s_{11} & s_{21} & \\ldots & s_{n_k1} \\\\\n",
    "s_{12} & s_{22} & \\cdots & s_{n_k2}\\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots \\\\\n",
    "s_{1n_k} & s_{2n_k} & \\cdots & s_{n_kn_k}\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf G = \n",
    "\\begin{pmatrix}\n",
    "P_0 m_1 & P_0m_2 & \\ldots & P_0 m_{n_k} \\\\\n",
    "0 & (\\nabla m_2, \\nabla m_2)_{0,E} & \\cdots & (\\nabla m_2, \\nabla m_{n_k})_{0,E}\\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots \\\\\n",
    "0 & (\\nabla m_{n_k}, \\nabla m_2)_{0,E} & \\cdots & (\\nabla m_{n_k}, \\nabla m_{n_k})_{0,E}\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathbf B =& \n",
    "\\begin{pmatrix}\n",
    "P_0\\varphi_1 & \\cdots & P_0\\varphi_{N_k}\\\\\n",
    "(\\nabla m_2, \\nabla\\varphi_1)_{0, E} & \\cdots & (\\nabla m_2, \\nabla\\varphi_{N_k})_{0, E}\\\\\n",
    "\\cdots & \\cdots & \\cdots \\\\\n",
    "(\\nabla m_{n_k}, \\nabla\\varphi_1)_{0, E} & \\cdots & (\\nabla m_{n_k}, \\nabla\\varphi_{N_k})_{0, E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "=& \\begin{pmatrix}\n",
    "P_0\\varphi_1 & \\cdots & P_0\\varphi_{N_k}\\\\\n",
    "-(\\Delta m_2,\\varphi_1)_{0, E} & \\cdots & -(\\Delta m_2,\\varphi_{N_k})_{0, E}\\\\\n",
    "\\cdots & \\cdots & \\cdots \\\\\n",
    "-(\\Delta m_{n_k}, \\varphi_1)_{0, E} & \\cdots & -(\\Delta m_{n_k},\\varphi_{N_k})_{0, E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "& +\\begin{pmatrix}\n",
    "0 & \\cdots & 0\\\\\n",
    "(\\frac{ \\partial m_2}{\\partial \\mathbf n},\\varphi_1)_{\\partial E} & \\cdots & (\\frac{\\partial m_2}{\\partial \\mathbf n},\\varphi_{N_k})_{\\partial E}\\\\\n",
    "\\cdots & \\cdots & \\cdots \\\\\n",
    "(\\frac{\\partial m_{n_k}}{\\partial \\mathbf n}, \\varphi_1)_{\\partial E} & \\cdots & (\\frac{\\partial m_{n_k}}{\\partial \\mathbf n}, \\varphi_{N_k})_{\\partial E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因为单项式空间 $\\mathcal M_k(E) \\subset \\mathcal V_k(E)$, 所以 $\\mathcal M_k(E)$ 中所有单项式 $m_\\alpha$ 可以用 $\\mathcal V_k(E)$ 的基函数表示出来, \n",
    "\n",
    "$$\n",
    "m_\\alpha = \\sum_{i=1}^{N_k} dof_i(m_\\alpha)\\varphi_i\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf D = \\begin{pmatrix}\n",
    "dof_1(m_1) & dof_1(m_2) & \\cdots & dof_1(m_{n_k}) \\\\\n",
    "dof_2(m_1) & dof_2(m_2) & \\cdots & dof_2(m_{n_k}) \\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots \\\\\n",
    "dof_{N_K}(m_1) & dof_{N_K}(m_2) & \\cdots & dof_{N_K}(m_{n_k}) \n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\Pi^\\nabla: \\mathcal V_k(E)\\rightarrow \\mathcal{P}_k(E)$ 的矩阵表示如下: \n",
    "\n",
    "$$\\mathbf \\Pi^\\nabla = \\mathbf G^{-1}\\mathbf B$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\Pi^{\\nabla} (\\varphi_1, \\varphi_2, \\cdots, \\varphi_{N_k}) =& (m_1, m_2, \\cdots, m_{n_k}) \\mathbf G^{-1}\\mathbf B\\\\\n",
    "= &(\\varphi_1, \\varphi_2, \\cdots, \\varphi_{N_K})\\mathbf D\\mathbf G^{-1}\\mathbf B\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "所以 $\\Pi^{\\nabla}_*: \\mathcal V_k(E)\\rightarrow \\mathcal V_k(E) $ 的矩阵表示如下:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf{\\Pi}^{\\nabla}_*=\\mathbf D\\mathbf G^{-1}\\mathbf B\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另外, 易证 $\\mathbf G = \\mathbf B\\mathbf D$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 单元刚度矩阵构造 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义在多边形 $E$ 上的 Laplace 算子的单元刚度矩阵，i.e.,\n",
    "$$\n",
    "(\\mathbf{K}_E)_{ij} = (\\nabla \\varphi_i, \\nabla \\varphi_j) \\quad i,j=1,...,N_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用投影算子$\\Pi^{\\nabla}$, 作分解：\n",
    "$$ \n",
    "\\varphi_i = \\Pi^{\\nabla} \\varphi_i + (I-\\Pi^{\\nabla})\\varphi_i\\\\\n",
    "\\varphi_j = \\Pi^{\\nabla} \\varphi_j + (I-\\Pi^{\\nabla})\\varphi_j\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "(\\mathbf K_{E})_{i j} =& (\\nabla \\Pi^{\\nabla} \\varphi_i, \\nabla \\Pi^{\\nabla} \\varphi_j) + (\\nabla(I-\\Pi^{\\nabla})\\varphi_i,\\nabla(I-\\Pi^{\\nabla})\\varphi_j)\\\\\n",
    "& + (\\nabla \\Pi^{\\nabla} \\varphi_i,\\nabla(I-\\Pi^{\\nabla})\\varphi_j)+(\\nabla(I-\\Pi^{\\nabla})\\varphi_i,\\nabla \\Pi^{\\nabla} \\varphi_j)\\\\\n",
    "= &(\\nabla \\Pi^{\\nabla} \\varphi_i, \\nabla \\Pi^{\\nabla} \\varphi_j) + (\\nabla(I-\\Pi^{\\nabla})\\varphi_i,\\nabla(I-\\Pi^{\\nabla})\\varphi_j)\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "这样就把虚单元法单元刚度矩阵的计算分解为可计算的两部分， 其中第一部分："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "&\\begin{pmatrix}\n",
    "(\\nabla\\Pi^\\nabla\\varphi_1, \\nabla\\Pi^\\nabla\\varphi_1) & (\\nabla\\Pi^\\nabla\\varphi_1, \\nabla\\Pi^\\nabla\\varphi_2) & \\cdots & (\\nabla\\Pi^\\nabla\\varphi_1, \\nabla\\Pi^\\nabla\\varphi_{N_k})\\\\\n",
    "(\\nabla\\Pi^\\nabla\\varphi_2, \\nabla\\Pi^\\nabla\\varphi_1) & (\\nabla\\Pi^\\nabla\\varphi_2, \\nabla\\Pi^\\nabla\\varphi_2) & \\cdots & (\\nabla\\Pi^\\nabla\\varphi_2, \\nabla\\Pi^\\nabla\\varphi_{N_k})\\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots \\\\\n",
    "(\\nabla\\Pi^\\nabla\\varphi_{N_k}, \\nabla\\Pi^\\nabla\\varphi_1) & (\\nabla\\Pi^\\nabla\\varphi_{N_k}, \\nabla\\Pi^\\nabla\\varphi_2) & \\cdots & (\\nabla\\Pi^\\nabla\\varphi_{N_k}, \\nabla\\Pi^\\nabla\\varphi_{N_k})\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "=&\\int_E \\nabla\\Pi^\\nabla \\begin{pmatrix} \\varphi_1\\\\ \\varphi_2 \\\\ \\vdots\\\\\\varphi_{n_k} \\end{pmatrix} \\left(\\nabla\\Pi^\\nabla \\begin{pmatrix} \\varphi_1\\\\ \\varphi_2 \\\\ \\vdots\\\\\\varphi_{n_k} \\end{pmatrix} \\right)^T\\mathrm d \\mathbf x \\\\\n",
    "=&\\int_E [\\mathbf \\Pi^\\nabla]^T \\begin{pmatrix}\n",
    "\\nabla m_1 \\\\ \\nabla m_2 \\\\\\vdots\\\\ \\nabla m_{n_k}\n",
    "\\end{pmatrix}\n",
    "(\\nabla m_1, \\nabla m_2, \\cdots,  \\nabla m_{n_k}) \\mathbf \\Pi^\\nabla \\mathrm d \\mathbf x \\\\\n",
    "= & [\\mathbf \\Pi^\\nabla]^T \\tilde{G} \\mathbf \\Pi^\\nabla\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf K_E^h = [\\mathbf \\Pi^\\nabla]^T \\tilde{G} \\mathbf \\Pi^\\nabla + (\\mathbf I - \\mathbf \\Pi^\\nabla_*)^T(\\mathbf I - \\mathbf \\Pi^\\nabla_*)\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $L^2$ 投影\n",
    "\n",
    "定义投影算子 $\\Pi^0:\\mathcal V_k(E)\\to \\mathcal P_k(E)$, 满足：\n",
    "\n",
    "$$\n",
    "(p_k, \\Pi^0 v_h) = (p_k, v_h), \\quad \\forall p_k \\in \\mathcal P_k(E)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\Pi^0 (\\varphi_1, \\varphi_2, \\cdots, \\varphi_{N_k}) =& (m_1, m_2, \\cdots, m_{n_k}) \\mathbf H^{-1}\\mathbf C\\\\\n",
    "= &(\\varphi_1, \\varphi_2, \\cdots, \\varphi_{N_K})\\mathbf D\\mathbf H^{-1}\\mathbf C\n",
    "\\end{aligned}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中 \n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf H = \\begin{pmatrix}\n",
    "(m_1, m_1) & (m_1, m_2) & \\cdots & (m_1, m_{n_k}) \\\\\n",
    "(m_2, m_1) & (m_2, m_2) & \\cdots & (m_2, m_{n_k}) \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\ \n",
    "(m_{n_k}, m_1) & (m_{n_k}, m_2) & \\cdots & (m_{n_k}, m_{n_k}) \n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf C = & \\begin{pmatrix}\n",
    "(m_1, \\varphi_1) & (m_1, \\varphi_2) & \\cdots & (m_1, \\varphi_{N_k}) \\\\\n",
    "(m_2, \\varphi_1) & (m_2, \\varphi_2) & \\cdots & (m_2, \\varphi_{N_k})\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
    "(m_{n_k}, \\varphi_1) & (m_{n_k}, \\varphi_2) & \\cdots & (m_{n_k}, \\varphi_{N_k})\n",
    "\\end{pmatrix}\\\\\n",
    "=& \\begin{pmatrix}\n",
    "\\mathbf C_{n_{k-2}\\times kN^V} & \\mathbf C_{n_{k-2}\\times n_{k-2}}\\\\\n",
    "\\mathbf C_{(n_{k} - n_{k-2})\\times kN^V} & \\mathbf C_{(n_{k} - n_{k-2})\\times n_{k-2}}\n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中:\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf C_{n_{k-2}\\times kN^V} = & \\mathbf 0 \\\\\n",
    "\\mathbf C_{n_{k-2}\\times n_{k-2}} = &|E| \\mathbf I\\\\\n",
    "(\\mathbf C_{(n_{k} - n_{k-2})\\times kN^V}, \\mathbf C_{(n_{k} - n_{k-2})\\times n_{k-2}}) = & \\mathbf H\\mathbf G^{-1}\\mathbf B[n_{k-2}+1:n_k,1:N_k]\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "C_{\\alpha i} = \n",
    "\\begin{cases}\n",
    "|E|, &i = kN^V+\\alpha\\\\\n",
    "(\\mathbf H\\mathbf G^{-1}\\mathbf B)_{\\alpha i}, &n_{k-2} < \\alpha \\leq n_k\\\\\n",
    "0, & \\text{others}\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记 $\\Pi^0: \\mathcal V_k(E) \\rightarrow \\mathcal{P}_k(E)$ 和 $\\Pi^0_*: \\mathcal V_k(E) \\rightarrow \\mathcal V_k(E)$:\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf \\Pi^0_* = \\mathbf D\\mathbf \\Pi^0 = \\mathbf D\\mathbf H^{-1} \\mathbf C \n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 质量矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "定义在多边形 $E$ 上的 Laplace 算子的质量矩阵，i.e.,\n",
    "$$\n",
    "(\\mathbf{M}_E)_{ij} = (\\varphi_i, \\varphi_j) \\quad i,j=1,...,N_k\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用投影算子$\\Pi^{\\nabla}$, 作分解：\n",
    "$$ \n",
    "\\varphi_i = \\Pi^{0} \\varphi_i + (I-\\Pi^{0})\\varphi_i\\\\\n",
    "\\varphi_j = \\Pi^{0} \\varphi_j + (I-\\Pi^{0})\\varphi_j\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "(\\mathbf M_{E})_{i j} =& (\\Pi^{0} \\varphi_i, \\Pi^{0} \\varphi_j) + ((I-\\Pi^{0})\\varphi_i,(I-\\Pi^{0})\\varphi_j)\\\\\n",
    "& + (\\Pi^{0} \\varphi_i,(I-\\Pi^{0})\\varphi_j)+((I-\\Pi^{0})\\varphi_i,\\Pi^{0} \\varphi_j)\\\\\n",
    "= &(\\Pi^{0} \\varphi_i,  \\Pi^{0} \\varphi_j) + ((I-\\Pi^{0})\\varphi_i,(I-\\Pi^{0})\\varphi_j)\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "这样就把虚单元法单元质量矩阵的计算分解为可计算的两部分， 其中第一部分："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "&\\begin{pmatrix}\n",
    "(\\Pi^0\\varphi_1, \\Pi^0\\varphi_1) & (\\Pi^0\\varphi_1, \\Pi^0\\varphi_2) & \\cdots & (\\Pi^0\\varphi_1, \\Pi^0\\varphi_{N_k})\\\\\n",
    "(\\Pi^0\\varphi_2, \\Pi^0\\varphi_1) & (\\Pi^0\\varphi_2, \\Pi^0\\varphi_2) & \\cdots & (\\Pi^0\\varphi_2, \\Pi^0\\varphi_{N_k})\\\\\n",
    "\\cdots & \\cdots & \\cdots & \\cdots \\\\\n",
    "(\\Pi^0\\varphi_{N_k}, \\Pi^0\\varphi_1) & (\\Pi^0\\varphi_{N_k}, \\Pi^0\\varphi_2) & \\cdots & (\\Pi^0\\varphi_{N_k}, \\Pi^0\\varphi_{N_k})\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "=&\\int_E \\Pi^0 \\begin{pmatrix} \\varphi_1\\\\ \\varphi_2 \\\\ \\vdots\\\\\\varphi_{n_k} \\end{pmatrix} \\left(\\Pi^0 \\begin{pmatrix} \\varphi_1\\\\ \\varphi_2 \\\\ \\vdots\\\\\\varphi_{n_k} \\end{pmatrix} \\right)^T\\mathrm d \\mathbf x \\\\\n",
    "=&\\int_E [\\mathbf \\Pi^0]^T \\begin{pmatrix}\n",
    " m_1 \\\\  m_2 \\\\\\vdots\\\\  m_{n_k}\n",
    "\\end{pmatrix}\n",
    "( m_1,  m_2, \\cdots,   m_{n_k}) \\mathbf \\Pi^0 \\mathrm d \\mathbf x \\\\\n",
    "= & [\\mathbf \\Pi^0]^T \\mathbf H \\mathbf \\Pi^0\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf M_E^h = \\mathbf C^T\\mathbf H^{-1}\\mathbf C + |E|(\\mathbf I - \\mathbf \\Pi^0_*)^T(\\mathbf I - \\mathbf \\Pi^0_*)\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 右端项\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf F = \\begin{pmatrix}\n",
    "(f, \\Pi^0_k\\varphi_1) \\\\\n",
    "(f, \\Pi^0_k\\varphi_2) \\\\\n",
    "\\vdots\\\\\n",
    "(f, \\Pi^0_k\\varphi_{N^{dof}})\n",
    "\\end{pmatrix}\n",
    "=(\\mathbf \\Pi^0)^T\\begin{pmatrix}\n",
    "(f, m_1)\\\\ (f, m_2) \\\\ \\vdots \\\\ (f, m_{n_k})\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 各种情形\n",
    "\n",
    "### $k= 1$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf G = &\n",
    "\\begin{pmatrix}\n",
    "P_0 m_1 & P_0m_2 &  P_0 m_{3} \\\\\n",
    "0 & (\\nabla m_2, \\nabla m_2)_{0,E} &  (\\nabla m_2, \\nabla m_{3})_{0,E}\\\\\n",
    "0 & (\\nabla m_{3}, \\nabla m_2)_{0,E} &  (\\nabla m_{3}, \\nabla m_{3})_{0,E}\n",
    "\\end{pmatrix}\\\\\n",
    "= & \\begin{pmatrix}\n",
    "1 & \\sum_{i=1}^{N^V} \\frac{x_i - x_{E}}{h_EN^V} & \\sum_{i=1}^{N^V} \\frac{y_i - y_{E}}{h_EN^V}\\\\\n",
    "0 & 1 & 0\\\\\n",
    "0 & 0 & 1\n",
    "\\end{pmatrix}\\\\\n",
    "=& \\begin{pmatrix}\n",
    "1 & 0 & 0\\\\\n",
    "0 & 1 & 0\\\\\n",
    "0 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf B =& \n",
    "\\begin{pmatrix}\n",
    "P_0\\varphi_1 & \\cdots & P_0\\varphi_{N^{V}}\\\\\n",
    "(\\nabla m_2, \\nabla\\varphi_1)_{0, E} & \\cdots & (\\nabla m_2, \\nabla\\varphi_{N^{V}})_{0, E}\\\\\n",
    "(\\nabla m_{3}, \\nabla\\varphi_1)_{0, E} & \\cdots & (\\nabla m_{3}, \\nabla\\varphi_{N^{V}})_{0, E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "=& \\begin{pmatrix}\n",
    "0 & \\cdots & 0\\\\\n",
    "-(\\Delta m_2,\\varphi_1)_{0, E} & \\cdots & -(\\Delta m_2,\\varphi_{N^{V}})_{0, E}\\\\\n",
    "-(\\Delta m_{3}, \\varphi_1)_{0, E} & \\cdots & -(\\Delta m_{3},\\varphi_{N^{V}})_{0, E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "& +\\begin{pmatrix}\n",
    "P_0\\varphi_1 & \\cdots & P_0\\varphi_{N^{V}}\\\\\n",
    "(\\frac{ \\partial m_2}{\\partial \\mathbf n},\\varphi_1)_{\\partial E} & \\cdots & (\\frac{\\partial m_2}{\\partial \\mathbf n},\\varphi_{N^{V}})_{\\partial E}\\\\\n",
    "(\\frac{\\partial m_{3}}{\\partial \\mathbf n}, \\varphi_1)_{\\partial E} & \\cdots & (\\frac{\\partial m_{3}}{\\partial \\mathbf n}, \\varphi_{N^{V}})_{\\partial E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "=&\\begin{pmatrix}\n",
    "\\frac{1}{N^V} & \\cdots & \\frac{1}{N^V}\\\\\n",
    "(\\frac{ \\partial m_2}{\\partial \\mathbf n},\\varphi_1)_{\\partial E} & \\cdots & (\\frac{\\partial m_2}{\\partial\\mathbf n},\\varphi_{N^{V}})_{\\partial E}\\\\\n",
    "(\\frac{\\partial m_{3}}{\\partial\\mathbf n}, \\varphi_1)_{\\partial E} & \\cdots & (\\frac{\\partial m_{3}}{\\partial\\mathbf n}, \\varphi_{N^{V}})_{\\partial E}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "=&\\begin{pmatrix}\n",
    "\\frac{1}{N^V} & \\cdots  & \\frac{1}{N^V}\\\\\n",
    "\\frac{1}{h_E}\\sum_e \\mathbf n_e\\int_e\\varphi_1& \\cdots& \\frac{1}{h_E}\\sum_e \\mathbf n_e\\int_e\\varphi_{N^V_P})\n",
    "\\end{pmatrix}\\\\\n",
    "=& \\begin{pmatrix}\n",
    "\\frac{1}{N^V} & \\cdots & \\frac{1}{N^V}\\\\\n",
    "\\frac{1}{h_E}\\mathbf d_1 & \\cdots & \\frac{1}{h_E} \\mathbf d_{N_E^V}\\\\\n",
    "\\end{pmatrix}\\\\\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "其中 $$\\mathbf d_{i} = \\frac{1}{2}\\mathbf W^T(\\mathbf x_{i+1} - \\mathbf x_{i-1})$$ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf D = &\\begin{pmatrix}\n",
    "dof_1(m_1) & dof_1(m_2) &  dof_1(m_{3}) \\\\\n",
    "dof_2(m_1) & dof_2(m_2) &  dof_2(m_{3}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "dof_{N^{v}}(m_1) & dof_{N^{v}}(m_2) & dof_{N^{v}}(m_{3}) \n",
    "\\end{pmatrix}\\\\\n",
    "=& \\begin{pmatrix}\n",
    "1 & \\frac{x_1 - x_{c,E}}{h_E} &  \\frac{y_1 - y_{c,E}}{h_E} \\\\\n",
    "1 & \\frac{x_2 - x_{c,E}}{h_E} &  \\frac{y_2 - y_{c,E}}{h_E} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "1 & \\frac{x_{N^V} - x_{c,E}}{h_E} & \\frac{y_{N^V} - y_{c,E}}{h_E}\n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\mathbf H = \\begin{pmatrix}\n",
    "(m_1, m_1) & (m_1, m_2) &  (m_1, m_{3}) \\\\\n",
    "(m_2, m_1) & (m_2, m_2) &  (m_2, m_{3}) \\\\\n",
    "(m_{3}, m_1) & (m_{3}, m_2) & (m_{3}, m_{3}) \n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\mathbf C =& \\begin{pmatrix}\n",
    "(m_1, \\varphi_1) & (m_1, \\varphi_2) & \\cdots & (m_1, \\varphi_{N^{V}}) \\\\\n",
    "(m_2, \\varphi_1) & (m_2, \\varphi_2) & \\cdots & (m_2, \\varphi_{N^{V}})\\\\\n",
    "(m_{3}, \\varphi_1) & (m_{3}, \\varphi_2) & \\cdots & (m_{3},\\varphi_{N^{V}})\n",
    "\\end{pmatrix}\\\\\n",
    "=& \\mathbf H\\mathbf G^{-1}\\mathbf B\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
